Deep dive regional communities

source('../env.R')

Species in communities

It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.

Load data

communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))

Rows: 2462 Columns: 9── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (1): city_id
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

communities_summary = communities %>% group_by(city_id) %>% summarise(
  regional_pool_size = n(), 
  urban_size_high = sum(present_urban_high), 
  urban_size_med = sum(present_urban_med), 
  urban_size_low = sum(present_urban_low)
)

communities_summary_long = bind_rows(
  communities_summary %>% rename(urban_pool_size = 'urban_size_high') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'High'),
  communities_summary %>% rename(urban_pool_size = 'urban_size_med') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'Medium'),
  communities_summary %>% rename(urban_pool_size = 'urban_size_low') %>% dplyr::select(city_id, regional_pool_size, urban_pool_size) %>% mutate(is_urban_threshold = 'Low')
)
communities_summary_long

city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))

Warning: st_centroid assumes attributes are constant over geometries of xWarning: st_centroid does not give correct centroids for longitude/latitude data

community_data_metrics = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics.csv')) %>%
  dplyr::select(city_id, mntd_normalised, fd_normalised, urban_pool_size, is_urban_threshold) %>%
  left_join(read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))) %>%
  left_join(communities_summary_long) %>%
  left_join(city_points[,c('city_id', 'city_nm')]) %>%
  rename(city_name='city_nm') %>%
  mutate(is_urban_threshold = factor(is_urban_threshold, levels = c('low', 'medium', 'high'), labels = c('Low', 'Medium', 'High'), ordered = T)) %>%
  arrange(city_id)

Rows: 792 Columns: 93── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): is_urban_threshold
dbl (92): city_id, mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fd_normalised, fd_actual, fd_min, fd_max, fd_mean, fd_sd, FRic_normal...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 342 Columns: 2── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id, urban_pool_size, is_urban_threshold)`Joining with `by = join_by(city_id)`

community_data_metrics

Load trait data

traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))

Rows: 304 Columns: 5── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (4): gape_width, trophic_trait, locomotory_trait, mass
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(traits)

fetch_normalised_traits = function(required_species_list) {
  required_traits = traits %>% filter(jetz_species_name %in% required_species_list)
  
  required_traits$gape_width_normalised = normalise(required_traits$gape_width, min(required_traits$gape_width), max(required_traits$gape_width))
  required_traits$trophic_trait_normalised = normalise(required_traits$trophic_trait, min(required_traits$trophic_trait), max(required_traits$trophic_trait))
  required_traits$locomotory_trait_normalised = normalise(required_traits$locomotory_trait, min(required_traits$locomotory_trait), max(required_traits$locomotory_trait))
  required_traits$mass_normalised = normalise(required_traits$mass, min(required_traits$mass), max(required_traits$mass))
  
  traits_normalised_long = required_traits %>% pivot_longer(cols = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), names_to = 'trait', values_to = 'normalised_value') %>% dplyr::select(jetz_species_name, trait, normalised_value)
  traits_normalised_long$trait = factor(traits_normalised_long$trait, levels = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), labels = c('Gape Width', 'Trophic Trait', 'Locomotory Trait', 'Mass'))
  
  traits_normalised_long
}

fetch_normalised_traits(c('Aplopelia_larvata', 'Chalcophaps_indica', 'Caloenas_nicobarica'))

Read in our phylogeny

phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')

Load resolve ecoregions

resolve = read_resolve()

Warning: st_buffer does not correctly buffer longitude/latitude datadist is assumed to be in decimal degrees (arc_degrees).
Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

Create helper functions

to_species_matrix = function(filtered_communities) {
  filtered_communities %>% 
    dplyr::select(city_id, jetz_species_name) %>% 
    distinct() %>%
    mutate(present = TRUE) %>% 
    pivot_wider(
      names_from = jetz_species_name, 
      values_from = "present", 
      values_fill = list(present = F)
    ) %>% 
    tibble::column_to_rownames(var='city_id')
}

community_nmds = function(filtered_communities) {
  species_matrix = to_species_matrix(filtered_communities)
  nmds <- metaMDS(species_matrix, k=2, trymax = 30)
  nmds_result = data.frame(scores(nmds))
  nmds_result$city_id = as.double(rownames(scores(nmds)))
  rownames(nmds_result) = NULL
  
  nmds_result
}

https://www.datacamp.com/tutorial/k-means-clustering-r

scree_plot = function(community_nmds_data) {
  # Decide how many clusters to look at
  n_clusters <- min(10, nrow(community_nmds_data) - 1)
  
  # Initialize total within sum of squares error: wss
  wss <- numeric(n_clusters)
  
  set.seed(123)
  
  # Look over 1 to n possible clusters
  for (i in 1:n_clusters) {
    # Fit the model: km.out
    km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
    # Save the within cluster sum of squares
    wss[i] <- km.out$tot.withinss
  }
  
  # Produce a scree plot
  wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
   
  scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
      geom_point(size = 4) +
      geom_line() +
      geom_hline(linetype="dashed", color = "orange", yintercept = wss) +
      scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
      xlab('Number of clusters')
  scree_plot
}

cluster_cities = function(city_nmds, cities_community_data, centers) {
  set.seed(123)
  kmeans_clusters <- kmeans(city_nmds[,c('NMDS1', 'NMDS2')], centers = centers, nstart = 20)
  city_nmds$cluster = kmeans_clusters$cluster
  cities_community_data %>% left_join(city_nmds) %>% mutate(cluster = as.factor(cluster))
}

plot_nmds_clusters = function(cluster_cities) {
  cluster_cities %>% dplyr::select(city_id, city_name, NMDS1, NMDS2, cluster) %>% distinct() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point() + geom_label_repel(aes(label = city_name))
}

plot_city_cluster = function(city_cluster_data_metrics, title) {
  species_in_cluster = communities %>% 
    filter(city_id %in% city_cluster_data_metrics$city_id) %>% 
    group_by(jetz_species_name, city_name, present_urban_high, present_urban_med, present_urban_low) %>% 
    summarise() %>%
    mutate(present_score = present_urban_high + present_urban_med + present_urban_low)

  species_in_cluster$present_score = factor(species_in_cluster$present_score, levels = c(0, 1, 2, 3), labels = c('Regional Pool', 'Low Threshold', 'Medium Threshold', 'High Threshold'))
  
  tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_in_cluster$jetz_species_name)))
    
  gg_tree = ggtree(tree_cropped)
  
  gg_presence = ggplot(species_in_cluster, aes(x=city_name, y=jetz_species_name)) + 
          geom_tile(aes(fill=present_score)) + scale_fill_manual(values = c("green", "yellow", "orange", "red")) + 
          theme_minimal() + xlab(NULL) + ylab(NULL) + 
          theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
          labs(fill='Presence')
  
  species_in_cluster_traits = fetch_normalised_traits(species_in_cluster$jetz_species_name)
  
  gg_traits = ggplot(species_in_cluster_traits, aes(x = trait, y = jetz_species_name, size = normalised_value)) + geom_point() + theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text.y=element_blank()) + xlab(NULL) + ylab(NULL) + labs(size = "Normalised Value")
  
  gg_cities_mntd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = mntd_normalised, fill = is_urban_threshold)) + geom_bar(stat = "identity", position = position_dodge(preserve = "single")) + scale_fill_manual(values = c("yellow", "orange", "red"))  + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("MNTD") + ylim(0, 1)
  
  gg_cities_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = fd_normalised, fill = is_urban_threshold)) + geom_bar(stat = "identity", position = position_dodge(preserve = "single")) + scale_fill_manual(values = c("yellow", "orange", "red")) + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("FD") + ylim(0, 1)
  
  gg_title = ggplot() + labs(title = title) + theme_minimal()
  
  gg_presence %>% insert_top(gg_cities_mntd, height = 0.5) %>% insert_top(gg_cities_fd, height = 0.5) %>% insert_left(gg_tree, width = 0.75) %>% insert_right(gg_traits, width = 0.5) %>% insert_top(gg_title, height = 0.06)
}

REGION_DEEP_DIVE_FIGURES_OUTPUT = mkdir(FIGURES_OUTPUT_DIR, 'appendix_regional_deep_dive')

Nearctic

nearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Nearctic')
nearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "San Jose"                 "Los Angeles"              "Concord"                  "Tijuana"                  "Bakersfield"             
 [6] "Fresno"                   "Sacramento"               "Mexicali"                 "Hermosillo"               "Las Vegas"               
[11] "Phoenix"                  "Tucson"                   "Durango"                  "Portland"                 "Chihuahua"               
[16] "Aguascalientes"           "Seattle"                  "Ciudad Juárez"            "San Luis Potosí"          "Mexico City"             
[21] "Saltillo"                 "Vancouver"                "Albuquerque"              "Monterrey"                "Nuevo Laredo"            
[26] "San Antonio"              "Austin"                   "Houston"                  "Dallas"                   "Oklahoma City"           
[31] "Calgary"                  "New Orleans"              "Kansas City"              "Omaha"                    "St. Louis"               
[36] "Bradenton"                "Tampa"                    "Minneapolis [Saint Paul]" "Atlanta"                  "Orlando"                 
[41] "Louisville"               "Chicago"                  "Indianapolis"             "Milwaukee"                "Virginia Beach"          
[46] "New York"                

nearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% nearctic_cities_community_data$city_id)) 

Run 0 stress 0.1161938 
Run 1 stress 0.1152814 
... New best solution
... Procrustes: rmse 0.03018951  max resid 0.1758497 
Run 2 stress 0.1152814 
... New best solution
... Procrustes: rmse 0.00004112701  max resid 0.0002340076 
... Similar to previous best
Run 3 stress 0.2542273 
Run 4 stress 0.1448163 
Run 5 stress 0.1198631 
Run 6 stress 0.1497522 
Run 7 stress 0.162113 
Run 8 stress 0.1152814 
... New best solution
... Procrustes: rmse 0.000008203033  max resid 0.00002624385 
... Similar to previous best
Run 9 stress 0.1145839 
... New best solution
... Procrustes: rmse 0.04424235  max resid 0.1777742 
Run 10 stress 0.1602924 
Run 11 stress 0.1148409 
... Procrustes: rmse 0.04440444  max resid 0.1744703 
Run 12 stress 0.1180524 
Run 13 stress 0.1176404 
Run 14 stress 0.1148408 
... Procrustes: rmse 0.04441218  max resid 0.1745631 
Run 15 stress 0.1145839 
... New best solution
... Procrustes: rmse 0.000007210488  max resid 0.00001868172 
... Similar to previous best
Run 16 stress 0.1180524 
Run 17 stress 0.1161938 
Run 18 stress 0.1148409 
... Procrustes: rmse 0.04440316  max resid 0.1744719 
Run 19 stress 0.114161 
... New best solution
... Procrustes: rmse 0.005102649  max resid 0.02607172 
Run 20 stress 0.1152814 
Run 21 stress 0.114161 
... Procrustes: rmse 0.000005376852  max resid 0.00001889389 
... Similar to previous best
*** Solution reached

nearctic_cities_nmds

scree_plot(nearctic_cities_nmds)

nearctic_cities = cluster_cities(city_nmds = nearctic_cities_nmds, cities_community_data = nearctic_cities_community_data, centers = 5)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(nearctic_cities)

nearctic_biomes = st_crop(resolve[resolve$REALM == 'Nearctic',c('REALM')], xmin = -220, ymin = 0, xmax = 0, ymax = 70)

although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries

 
ggplot() + 
  geom_sf(data = nearctic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = nearctic_cities, aes(geometry = geometry, color = cluster))

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_clusters.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 1

nearctic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Neartic cluster 1')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 2

nearctic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Neartic cluster 2')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 3

nearctic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Neartic cluster 3')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster3.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 4

nearctic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Neartic cluster 4')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 5

nearctic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Neartic cluster 5')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster5.jpg'))

Saving 7.29 x 4.51 in image

Neotropic

neotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Neotropic')
neotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "Culiacán"                  "Guadalajara"               "Morelia"                   "Acapulco"                  "Querétaro"                
 [6] "Cuernavaca"                "Puebla"                    "Oaxaca"                    "Xalapa"                    "Veracruz"                 
[11] "Tuxtla Gutiérrez"          "Villahermosa"              "Guatemala City"            "San Salvador"              "San Pedro Sula"           
[16] "Mérida"                    "Tegucigalpa"               "Managua"                   "San José"                  "Cancún"                   
[21] "Guayaquil"                 "Chiclayo"                  "Panama City"               "Trujillo"                  "Quito"                    
[26] "Havana"                    "Cali"                      "Lima"                      "Pereira"                   "Miami"                    
[31] "Medellín"                  "Ibagué"                    "Cartagena"                 "Kingston"                  "Bogota"                   
[36] "Barranquilla"              "Bucaramanga"               "Cúcuta"                    "Maracaibo"                 "Arequipa"                 
[41] "Barquisimeto"              "Santo Domingo"             "Maracay"                   "El Alto [La Paz]"          "Caracas"                  
[46] "Cochabamba"                "Viña del Mar [Valparaíso]" "Río Piedras [San Juan]"    "Barcelona"                 "Concepción"               
[51] "Santiago"                  "Mendoza"                   "Salta"                     "Cordoba"                   "Asuncion"                 
[56] "Buenos Aires"              "La Plata"                  "Ciudad del Este"           "Montevideo"                "Mar del Plata"            
[61] "Porto Alegre"              "São Paulo"                 "Santos"                    "Sao Jose dos Campos"      

neotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% neotropic_cities_community_data$city_id)) 

Run 0 stress 0.134619 
Run 1 stress 0.141222 
Run 2 stress 0.1348238 
... Procrustes: rmse 0.01072062  max resid 0.05476751 
Run 3 stress 0.1348237 
... Procrustes: rmse 0.01071747  max resid 0.05475478 
Run 4 stress 0.134433 
... New best solution
... Procrustes: rmse 0.006824883  max resid 0.04575085 
Run 5 stress 0.134619 
... Procrustes: rmse 0.006831176  max resid 0.04573137 
Run 6 stress 0.1405634 
Run 7 stress 0.1348046 
... Procrustes: rmse 0.006552356  max resid 0.04609573 
Run 8 stress 0.141222 
Run 9 stress 0.1346226 
... Procrustes: rmse 0.007200092  max resid 0.04572692 
Run 10 stress 0.1414222 
Run 11 stress 0.1363758 
Run 12 stress 0.1363675 
Run 13 stress 0.1344331 
... Procrustes: rmse 0.00002617473  max resid 0.00009028585 
... Similar to previous best
Run 14 stress 0.1344509 
... Procrustes: rmse 0.00118731  max resid 0.007067449 
... Similar to previous best
Run 15 stress 0.134619 
... Procrustes: rmse 0.006829458  max resid 0.04574559 
Run 16 stress 0.134433 
... Procrustes: rmse 0.00002418546  max resid 0.00008384028 
... Similar to previous best
Run 17 stress 0.134433 
... New best solution
... Procrustes: rmse 0.000008862921  max resid 0.00003722846 
... Similar to previous best
Run 18 stress 0.141222 
Run 19 stress 0.134433 
... Procrustes: rmse 0.00003797808  max resid 0.0001268688 
... Similar to previous best
Run 20 stress 0.1348046 
... Procrustes: rmse 0.006553689  max resid 0.04610243 
*** Solution reached

neotropic_cities_nmds

scree_plot(neotropic_cities_nmds)

neotropic_cities = cluster_cities(city_nmds = neotropic_cities_nmds, cities_community_data = neotropic_cities_community_data, centers = 5)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(neotropic_cities)

neotropic_biomes = resolve[resolve$REALM == 'Neotropic',c('REALM')]
 
ggplot() + 
  geom_sf(data = neotropic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = neotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_clusters.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 1

neotropic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Neotropic cluster 1')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 2

neotropic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Neotropic cluster 2')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 3

neotropic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Neotropic cluster 3')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster3.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 4

neotropic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Neotropic cluster 4')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 5

neotropic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Neotropic cluster 5')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster5.jpg'))

Saving 7.29 x 4.51 in image

Palearctic

palearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Palearctic')
palearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "Lisbon"                "Porto"                 "Marrakesh"             "Seville"               "Dublin"                "Málaga"               
 [7] "Madrid"                "Glasgow"               "Bilbao"                "Liverpool"             "Bristol"               "Manchester"           
[13] "Birmingham"            "Leeds"                 "Newcastle upon Tyne"   "Sheffield"             "Nottingham"            "Valencia"             
[19] "London"                "Toulouse"              "Paris"                 "Barcelona"             "Rotterdam [The Hague]" "Brussels"             
[25] "Amsterdam"             "Lyon"                  "Marseille"             "Dusseldorf"            "Nice"                  "Frankfurt am Main"    
[31] "Zurich"                "Oslo"                  "Stuttgart"             "Hamburg"               "Genoa"                 "Nuremberg"            
[37] "Copenhagen"            "Munich"                "Berlin"                "Dresden"               "Rome"                  "Prague"               
[43] "Stockholm"             "Poznan"                "Vienna"                "Wroclaw"               "Zagreb"                "Gdansk"               
[49] "Budapest"              "Krakow"                "Warsaw"                "Helsinki"              "Riga"                  "Belgrade"             
[55] "Lviv"                  "Sofia"                 "Thessaloniki"          "Minsk"                 "Athens"                "Kyiv"                 
[61] "Istanbul"              "Odesa"                 "Samsun"                "Luxor"                 "Tel Aviv"              "Jerusalem"            
[67] "Tbilisi"               "Yerevan"               "Abu Dhabi"             "Dubai"                 "Bishkek"              

palearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% palearctic_cities_community_data$city_id)) 

Run 0 stress 0.04203818 
Run 1 stress 0.05064867 
Run 2 stress 0.07206189 
Run 3 stress 0.04240445 
... Procrustes: rmse 0.03948386  max resid 0.2064496 
Run 4 stress 0.1164521 
Run 5 stress 0.05197845 
Run 6 stress 0.1332585 
Run 7 stress 0.05064862 
Run 8 stress 0.08350482 
Run 9 stress 0.04984194 
Run 10 stress 0.065365 
Run 11 stress 0.07340577 
Run 12 stress 0.04203776 
... New best solution
... Procrustes: rmse 0.00009496175  max resid 0.0003265871 
... Similar to previous best
Run 13 stress 0.05819636 
Run 14 stress 0.05530089 
Run 15 stress 0.1107012 
Run 16 stress 0.04223506 
... Procrustes: rmse 0.04136181  max resid 0.2169081 
Run 17 stress 0.1026199 
Run 18 stress 0.04223585 
... Procrustes: rmse 0.04135834  max resid 0.2172361 
Run 19 stress 0.08634193 
Run 20 stress 0.07085973 
*** Solution reached

palearctic_cities_nmds

scree_plot(palearctic_cities_nmds)

Warning: did not converge in 10 iterationsWarning: did not converge in 10 iterationsWarning: did not converge in 10 iterationsWarning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3550)Warning: did not converge in 10 iterationsWarning: did not converge in 10 iterationsWarning: Quick-TRANSfer stage steps exceeded maximum (= 3550)

palearctic_cities = cluster_cities(city_nmds = palearctic_cities_nmds, cities_community_data = palearctic_cities_community_data, centers = 6)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(palearctic_cities)

palearctic_biomes = resolve[resolve$REALM == 'Palearctic',c('REALM')]
 
ggplot() + 
  geom_sf(data = palearctic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = palearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_clusters.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 1

palearctic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Palearctic cluster 1')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 2

palearctic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Palearctic cluster 2')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 3

palearctic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Palearctic cluster 3')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster3.jpg'))

Saving 14 x 6 in image

Palearctic Cluster 4

palearctic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Palearctic cluster 4')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 5

palearctic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Palearctic cluster 5')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 6

palearctic_cities %>% filter(cluster == 6) %>% plot_city_cluster('Palearctic cluster 6')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster6.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic

afrotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Afrotropic')
afrotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
[1] "Cape Town"    "Johannesburg" "Pretoria"     "Kigali"       "Kampala"      "Arusha"       "Nairobi"      "Addis Ababa" 

afrotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% afrotropic_cities_community_data$city_id)) 

Run 0 stress 0.06095916 
Run 1 stress 0.06095916 
... New best solution
... Procrustes: rmse 0.00001059092  max resid 0.00001521966 
... Similar to previous best
Run 2 stress 0.1254289 
Run 3 stress 0.126785 
Run 4 stress 0.126785 
Run 5 stress 0.06095935 
... Procrustes: rmse 0.0001624829  max resid 0.0003126969 
... Similar to previous best
Run 6 stress 0.3007162 
Run 7 stress 0.126785 
Run 8 stress 0.126785 
Run 9 stress 0.06095923 
... Procrustes: rmse 0.0006753449  max resid 0.001228519 
... Similar to previous best
Run 10 stress 0.06095916 
... Procrustes: rmse 0.000006448791  max resid 0.00001242034 
... Similar to previous best
Run 11 stress 0.06095923 
... Procrustes: rmse 0.00004902177  max resid 0.00008420442 
... Similar to previous best
Run 12 stress 0.1246142 
Run 13 stress 0.126785 
Run 14 stress 0.126785 
Run 15 stress 0.1973021 
Run 16 stress 0.1246142 
Run 17 stress 0.06095904 
... New best solution
... Procrustes: rmse 0.0002457047  max resid 0.0004426671 
... Similar to previous best
Run 18 stress 0.0609592 
... Procrustes: rmse 0.0002974389  max resid 0.0005449775 
... Similar to previous best
Run 19 stress 0.1354466 
Run 20 stress 0.126785 
*** Solution reached

afrotropic_cities_nmds

scree_plot(afrotropic_cities_nmds)

afrotropic_cities = cluster_cities(city_nmds = afrotropic_cities_nmds, cities_community_data = afrotropic_cities_community_data, centers = 4)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(afrotropic_cities)

afrotropic_biomes = resolve[resolve$REALM == 'Afrotropic',c('REALM')]
 
ggplot() + 
  geom_sf(data = afrotropic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = afrotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_clusters.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic Cluster 1

afrotropic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Afrotropic cluster 1')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic Cluster 2

afrotropic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Afrotropic cluster 2')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic Cluster 3

afrotropic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Afrotropic cluster 3')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster3.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic Cluster 4

afrotropic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Afrotropic cluster 4')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan

indomalayan_cities_community_data = community_data_metrics %>% filter(core_realm == 'Indomalayan')
indomalayan_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
  [1] "Swat City"           "Srinagar"            "Jamnagar"            "Jammu"               "Rajkot"              "Bikaner"             "Jodhpur"            
  [8] "Jalandhar"           "Ahmedabad"           "Bhavnagar"           "Ludhiana"            "Anand"               "Udaipur"             "Surat"              
 [15] "Vadodara"            "Ajmer"               "Chandigarh"          "Vasai-Virar"         "Mumbai"              "Jaipur"              "Delhi [New Delhi]"  
 [22] "Nashik"              "Dehradun"            "Kota"                "Pune"                "Haridwar"            "Dhule"               "Ujjain"             
 [29] "Indore"              "Ahmadnagar"          "Kolhapur"            "Jalgaon"             "Agra"                "Aurangabad"          "Sangli"             
 [36] "Belagavi"            "Gwalior"             "Budaun"              "Bareilly"            "Dharwad"             "Bhopal"              "Bhind"              
 [43] "Mangaluru"           "Solapur"             "Vijayapura"          "Akola"               "Latur"               "Kannur"              "Davanagere"         
 [50] "Thalassery"          "Amravati"            "Kalaburagi"          "Kozhikode"           "Guruvayur"           "Malappuram"          "Lucknow"            
 [57] "Thrissur"            "Mysuru"              "Kochi"               "Nagpur"              "Kollam"              "Jabalpur"            "Ettumanoor"         
 [64] "Hyderabad"           "Coimbatore"          "Bengaluru"           "Thiruvananthapuram"  "Tiruppur"            "Faizabad"            "Erode"              
 [71] "Prayagraj"           "Pratapgarh"          "Salem"               "Dindigul"            "Madurai"             "Tiruchirappalli"     "Durg"               
 [78] "Vellore"             "Tirupati"            "Raipur"              "Bilaspur"            "Vijayawada"          "Puducherry"          "Chennai"            
 [85] "Kathmandu"           "Colombo"             "Rajamahendravaram"   "Patna"               "Kandy"               "Bihar Sharif"        "Visakhapatnam"      
 [92] "Ranchi"              "Brahmapur"           "Jamshedpur"          "Darjeeling"          "Siliguri"            "Cuttack"             "Bhubaneshwar"       
 [99] "Jalpaiguri"          "Berhampore"          "Kolkata"             "Krishnanagar"        "Guwahati [Dispur]"   "Agartala"            "Silchar"            
[106] "Dimapur"             "Bangkok"             "George Town"         "Kuala Lumpur"        "Phnom Penh"          "Singapore"           "Hong Kong"          
[113] "Sha Tin"             "Hsinchu"             "Taichung"            "New Taipei [Taipei]" "Tainan"              "Denpasar"            "Kaohsiung"          
[120] "Kota Kinabalu"      

indomalayan_cities_nmds = community_nmds(communities %>% filter(city_id %in% indomalayan_cities_community_data$city_id)) 

Run 0 stress 0.1199339 
Run 1 stress 0.1221765 
Run 2 stress 0.1574272 
Run 3 stress 0.165272 
Run 4 stress 0.1180044 
... New best solution
... Procrustes: rmse 0.03780687  max resid 0.2489361 
Run 5 stress 0.1223128 
Run 6 stress 0.128065 
Run 7 stress 0.1462967 
Run 8 stress 0.1177488 
... New best solution
... Procrustes: rmse 0.01878365  max resid 0.1301728 
Run 9 stress 0.1480994 
Run 10 stress 0.1205201 
Run 11 stress 0.1380584 
Run 12 stress 0.1377785 
Run 13 stress 0.1241019 
Run 14 stress 0.1174802 
... New best solution
... Procrustes: rmse 0.003927371  max resid 0.03013146 
Run 15 stress 0.1617252 
Run 16 stress 0.1645131 
Run 17 stress 0.1603814 
Run 18 stress 0.1241049 
Run 19 stress 0.1596718 
Run 20 stress 0.1228515 
Run 21 stress 0.1597143 
Run 22 stress 0.1547365 
Run 23 stress 0.1488411 
Run 24 stress 0.1464657 
Run 25 stress 0.1441819 
Run 26 stress 0.1221274 
Run 27 stress 0.1195283 
Run 28 stress 0.115795 
... New best solution
... Procrustes: rmse 0.02711525  max resid 0.2451229 
Run 29 stress 0.1224208 
Run 30 stress 0.1167315 
*** No convergence -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    27: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin

indomalayan_cities_nmds

scree_plot(indomalayan_cities_nmds)

indomalayan_cities = cluster_cities(city_nmds = indomalayan_cities_nmds, cities_community_data = indomalayan_cities_community_data, centers = 6)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(indomalayan_cities)

indomalayan_biomes = resolve[resolve$REALM == 'Indomalayan',c('REALM')]
 
ggplot() + 
  geom_sf(data = indomalayan_biomes, aes(geometry = geometry)) + 
  geom_sf(data = indomalayan_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_clusters.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 1

indomalayan_cities %>% filter(cluster == 1) %>% plot_city_cluster('Indomalayan cluster 1')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 2

indomalayan_cities %>% filter(cluster == 2) %>% plot_city_cluster('Indomalayan cluster 2')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 3

indomalayan_cities %>% filter(cluster == 3) %>% plot_city_cluster('Indomalayan cluster 3')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster3.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 4

indomalayan_cities %>% filter(cluster == 4) %>% plot_city_cluster('Indomalayan cluster 4')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 5

indomalayan_cities %>% filter(cluster == 5) %>% plot_city_cluster('Indomalayan cluster 5')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster5.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 6

indomalayan_cities %>% filter(cluster == 6) %>% plot_city_cluster('Indomalayan cluster 6')

`summarise()` has grouped output by 'jetz_species_name', 'city_name', 'present_urban_high', 'present_urban_med'. You can override using the `.groups` argument.

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster6.jpg'))

Saving 7.29 x 4.51 in image